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Abstract 

We develop and analyse a discrete, one-dimensional model of cell motility which incorporates 
the effects of volume filling, cell-to-cell adhesion and chemotaxis. The formal continuum limit 
of the model is a nonlinear generalisation of the parabolic-elliptic Keller-Segel equations, with 
a diffusivity which can become negative if the adhesion coefficient is large. The consequent 
ill-posedness results in the appearance of spatial oscillations and the development of plateaus in 
numerical solutions of the underlying discrete model. A global-existence result is obtained for 
the continuum equations in the case of favourable parameter values and data, and a steady-state 
analysis which, amongst other things, accounts for high-adhesion plateaus is carried out. For 
ill-posed cases, a singular Stefan-problem formulation of the continuum limit is written down 
and solved numerically, and the numerical solutions are compared with those of the original 
discrete model. 

1 Introduction 

Of late, there has been considerable interest in formulating continuum models for cell structures 
generated by cell-to-cell adhesion, the motivation being to facilitate mathematical analysis and 
efficient numerical simulation of processes such as de novo blood-vessel synthesis (i.e. vasculogenesis) 
and cancer invasion, for example. Two recent attempts in this direction were made, respectively, by 
Armstrong et al. [3], who began with a nonlocal integro-differential equation in which the kernel is 
integrated over a given cell-sensing radius, and by Anguige & Schmeiser J7, who wrote down a simple 
random-walk model accounting for adhesion, diffusion and volume filling. For both approaches, it 
turns out that the limiting macroscopic model (a nonlinear parabolic equation) can be ill posed if 
the adhesion is sufficiently strong, which leads to interesting pattern-forming behaviour in solutions 
of the underlying microscopic models, but which also makes mathematical analysis rather more 
difficult than one would like [2 [5] . 

Our intention in this paper is to extend the modelling and analysis of [U [2] in a rather obvious 
way, namely, by factoring in the directed response of cells to an extracellular chemical gradient (i.e. 
chemotaxis), and then examining the resulting interaction between such long-range signalling and 
the short-range signalling of cell-to-cell adhesion. 

* e-mail: keith. anguigeCSoeaw. ac. at 
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Following the approach previously adopted in [5] , our 1-d model for cell adhesion, diffusion and 
chemotaxis will take the form of the random-walk system 

^ = 7-lip._i + %-,p,+^ - [%+ + T-)P^ , (1) 

for the approximate cell densities pi G [0, 1], on a uniform grid of points Xi — ih, the quantities 7^^ 
being the transitional probabilities per unit time of a one-step jump from i to i ± 1. 

Taking inspiration from [2], as well as [7j, the scaled transitional probabilities are chosen to be 

^^ = ^(l-mi)(l-«ATi)(l + y(^^±i-^0) , (2) 

where a G [0, 1] is the adhesion coefficient, xo G [0, oo) the chemotactic sensitivity, and Si the 
concentration of chemoattractant at the point Xi. Here, the first factor in parentheses models 
volume filling, the second models adhesion, and the third chemotaxis. 

The expression ([2]) includes, as special cases, both our previous model for adhesion/diffusion 
(Xo = 0) [2], and the model for linear diffusion, chemotaxis and volume filling (a = 0) presented in 
0. 

Upon taking the continuum limit of (Il])-([2]), by means of Taylor expansion in powers of h, one 
ends up with the advection-diffusion equation 

9p d f dp dS\ 



dt dx \ dx dx 
where the diffusivity is given by 

/ 2\ ^ 4 
D{p)=ia[^p--J +l-3«, (4) 

and the chemotactic-sensitivity function by 

x(p)=Xo(l-p)(l-«p). (5) 

The continuum equation for the chemoattractant is taken to be the usual non-dimensionalised 
quasi-steady-state equation 

AS' = 5 - p. (6) 

Equations ©-([G]) constitute a nonlinear generalisation of the classical Keller-Segel chemotaxis 
model, depending on just two parameters, xo and a. They are to be solved on the domain — (0, L), 
subject to the Neumann boundary conditions on both p and S. 

Note that the presence of the volume-filling term in x{p) entails that ([3]) has the following 
Maximum Principle: 

< /5(0) < 1 ^ < p{t) < 1. (7) 
Moreover, ([6]) implies a Maximum Principle for S: 

< p{t) < 1 ^ < S{t) < 1. (8) 

As in [2], however, ^ can be ill-posed if a > |, since in that case there is an interval of values 
of p, centred around p = |, for which the diffusivity is negative. Explicitly, we have from (|4]) that 
D{p) < whenever 



pel^:=l , = (p , p«), (9) 
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and D{p) > otherwise. Note that the width of la increases a.s a y 1, and that /i — (i, l). 

Equation ([3]), therefore, is certainly ill posed if the initial density profile hits /„, and we saw in 
[5] that this ill-posedness is related to the presence of oscillations and the formation of plateaus in 
solutions of ID) in the special case xo = 0. As a consequence of chemotactic aggregation, however, 
we might also expect (IS])-® to be ill posed for small initial densities, provided xo is large enough 
and a > |. 

What interests us most in this paper, then, is possible singular pattern-forming behaviour in 
the case a > f and xo > 0. This turns out to be rather difficult to investigate analytically, but 
considerable insight may nevertheless be gained with the aid of numerical simulations. 

The paper is organised as follows. 

In Section 2 we carry out a steady-state analysis of (l3])-([6|), and this is followed in Section 3 by 
some global-existence results for favourable parameter values and initial data. 

In Section 4 we report on numerical simulations of ©-([I]), which show that singular (i.e. sharp- 
edged) aggregation patterns can be generated by small data, provided the chemotactic sensitivity 
and the adhesion coefficient are chosen large enough. 

Finally, and in analogy with [1], we consider the idea that a Stefan-problem- type framework, in 
which the density is allowed to jump across the unstable region la, might be an appropriate way of 
treating ©-([I]) as the limit of ^ in ill-posed cases. Although such problems seem to be analytically 
intractable at present, there is a sense in which solving them numerically may nevertheless be more 
efficient than discretising the Neumann problem for ((21)-® directly, since this (the direct approach) 
requires a very fine mesh to properly resolve the singular behaviour observed near la', simulations 
obtained using both of these methods are compared and contrasted in Section 5. 



2 Steady-state analysis 

We now show that essentially the same techniques employed in the case of linear diffusion (see, e.g., 
[7]) can also be used to investigate steady states of ([HI)-®. For a > |, it is also possible, as in the 
case Xo = p], to construct discontinuous weak solutions in which p has finitely many jumps across 
the mistablc region la- 



2.1 Linear stability of uniform steady states 

Linearising ®, ®, ®, ® around a uniform steady state p = S = p, and performing a standard 
stability calculation, gives the dispersion relation 

A = ^ ^\TJ , (10) 



L2 \^ ^'^^ 2.2+^7^2^ 

for the growth rate A and the wave number k. 

Thus, if Xo is so small that x{p)p < D{p), then A < V/c, and the uniform steady state is 
unconditionally stable. If, on the other hand, x{p)p > D{p), which necessarily occurs when a > | 
and p d la, for example, then A > for small wave numbers, and the uniform steady state is thus 
unstable to long-wavelength perturbations. 

It is also elementary to show that the dominant wavemode is determined by 



kT^Y X{p)p 



L J V D{p) 



(11) 
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2.2 Global L°°-stability of uniform steady states for a < | and xo small 

For a given mass M, Q-® always has the uniform solution (p, S) — {j5,j5), where p = M/L. One 
expects such a solution to be a global attractor provided a < | and xo is sufficiently small, and this 
is the content of the following theorem. 

Theorem 1 // 

Xominfl,yi72} niax (l-p)(l-«p)p< (12) 

^. J 0<p<l o 

then any smooth, global solution pair {p, S) of ^-(d^ satisfies \\p ~ pWooit) < ci — p||2(0)e~'^^* and 
\\S — p\\H'^{t) < \\p ~ p\\2{0)e~'^'^'^ , for some positive constants Ci, C2- 

Proof. First of all, in we subtract p from p, multiply through by p — p, and integrate by 
parts to get 

"P-P\\2 = - D{p){{p- p)xf dx + xo 1 {p-p)x{^-p){^-ap)pSj,dx. (13) 



2dt 

Next, note that squaring ([6|) gives, with the aid of an integration by parts 



{S-pY + 2{{S-p)xY + {{S~p)xxY dx^ {p-pYdx, (14) 

Jo 

and that differentiating ([6|) and carrying out the same procedure gives 

{{S - p)xY + 2((^ - pYx? + {{S - pY.x f dx^ [ {{p - pY? dx, (15) 



since the Neumann conditions kill the boundary terms in both cases. 

In particular, and (jlSp. together with the Poincare inequality, imply that 



||(^ - pYh < min {l, /l72} \\{P - pYh- (16) 
Hence, condition (|12p and a further application of the Poincare inequality imply that 

Ij^\\p-P\\l < -e|l(P"P)x||2 

< -^\\P-P\\l (17) 

for some e > 0, and therefore also 

||p-p||2(^)<llp-p||2(0)e-*/^ (18) 

and, by 

WS^pWh^ < ||p-p|l2(0)e-^*/^. (19) 
Next, for a given n G N, we integrate the first of pT|) from t = n to n + 1, thus obtaining 

n+l 

||(p-p).||^di<C||p-p||2(0)e— /^ (20) 



/ 
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which impHes that Vn e N, 3 C" e [n,n+ 1] such that ||(p - p)xhiC) < C\\p - p||2(0)e-<="/-^. 
Therefore, by Sobolev imbedding and L^-decay (or Poincare again), 

\\p-p\UC)<C\\p-p\me-^^-/^. (21) 

From here, we proceed with a comparison argument, which will demonstrate that p — p can grow 
(pointwise) at worst exponentially between the times C", which of course satisfy |C" — C"^^l ^ 2. 
Expanding ([3]), and using ^ to substitute for the Laplacian, we see that u :— p — p satisfies 

du d f ^, .du\ d , . dudS , , , „ xn 

Tt - (^(" + ^) j - d-p^^'^^'^d-.^. ^ + - 

Thus, by (|19p . we have that u is a subsolution of the following problem, for t £ [C", C"^^)- 

^ " ^ (^^^' + + ^^'^ ^"^ " ^ll2(0)e"'*/^ (23) 

subject to 

wMt)^w^{L,t)^Q, w;(x,C) = CI|p-p||2(0)e-^^"/^, (24) 

where B{x,t) is a bounded fmiction, say \B{x,t)\ < Co V x,t, and A(x,t), D{x,t) are smooth 
functions, with D{x,t) > 0. 

It is easy to check that w := Ce""^"/^ ||p-p||2(0)e'^i(*~^") is a supersolution of provided 
Ci > Co is chosen sufficiently large. 

Subsolutions can be constructed similarly, and we therefore arrive at 

||^^||oo(i)<C||p-p||2(0)e-^*/^, (25) 

Vt > 0, as required □ 

Clearly the argument used in the proof of Theorem 1 can be bootstrapped to obtain 

Theorem 2 The uniform steady state. {p,S) = {p,S), of {3^-® is a local L°° -attractor provided 

min (l, /l72) X{P)P < D{p). (26) 

Note that condition (pS)) is consistent with the linear stability analysis of Section 2.1, and that 
this condition holds for xo sufficiently small, provided either a < | or both a > | and p ^ /„. 

2.3 Non-trivial steady states; dynamical systems 
2.3.1 The case a < | 

We are looking for smooth solutions of 

Sxx ~ S + p = 0, 

subject to the boundary conditions Sx = Px = x — 0, L. 
First of all, equation ([27]) can be integrated to get 

G{p) ^S^C, 
5 



(27) 
(28) 

(29) 



where G is a primitive for D{p)/{x{p)p) {G strictly increasing since a < |), and C is a constant of 
integration. 

Using (pS)) . this gives us 

G{p)^^^G{p)~ p + C, (30) 

or, in terms of cr :— G{p), 

(T.. =a-G-i(a) + C, (31) 

which is a Hamiltonian dynamical system, whose critical points are necessarily saddles or centres. 

It is easy to see that pip has either one or three critical points, depending on a,xo and the 
choice of C. If xo is small, then there is just one critical point, and the only possible steady states 
satisfying the boundary conditions are uniform solutions. If, on the other hand, xo is sufficiently 
large, then the function G{p) ~ p has two extrema, and for an interval of values of C there are three 
critical points of (j3ip . comprising two saddles and a nonlinear centre inbetween. Solutions of the 
Neumann problem for (|3ip are then realised as half-orbits around the centre, or integer multiples 
thereof. 

If we have obtained a solution, p, of ((30| , and if S is then determined by the Neumann problem 
for ((28|) . then it is easy to see that the original steady-state equation ([27|) also holds. Indeed, using 
([28]) to substitute for the linear p-term in the rhs of ((30|) leads to the energy equality 

\\G{p)-S + C\\m=Q, (32) 

which gives the desired result, since G'{p) = D{p)/{x{p)p)- 

By inspecting G{p), it is also easy to see that there must be sufficient (and also not too much) 
mass for such solutions to exist. 

In addition, for a given C there is a minimum value of the domain length which allows for 
non- trivial steady-state solutions, and this is given by 



where Pc is the location of the nonlinear centre; L* is calculated by linearising about the centre. 

In fact, there is a critical curve in (a,xo)-space which divides the parameter region in which 
there is the possibility of three critical points from that in which there can only be a single one. It 
is obtained by solving F{p) = F'{p) = 0, where F{p) := D{p) — x{p)P: ^i^d is plotted in Figure 1 as 
the boundary between regions (ii) and (iii). 

A fuller explanation of Figure 1 is as follows. Region (i) is where there are no non-uniform steady 
states, and where uniform steady states have been proved to be nonlinearly stable, no matter how 
large L is (see Theorem 1). In region (ii) there are still no non-uniform steady states, and uniform 
steady states are always linearly stable; nonlinear stability has only been proved for small enough 
L (Theorem 1 again). In region (iii), non- uniform steady states become possible, for appropriate G 
and large enough L, and uniform states can lose their stability. In region (iv), the diffusivity can 
turn negative, and one may look for both smooth and discontinuous steady-state solutions, as we 
now discuss. 

2.3.2 The case a > |; weak solutions 

For a > |, non-uniform smooth solutions which miss are possible for all choices of xo > Oi 
provided L is sufficiently large. For small xo these solutions have to be small-amplitude oscillations 
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Figure 1: (a, xo)-space. 
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just outside the unstable region /„, while for large xo they have to be small-amplitude oscillations 
near p = or p = 1. These conclusions are again reached by considering the form of G{p). 

We can also look for weak solutions {p, S) of P7)) - (P5|) such that p contains finitely many jumps 
across la, but is smooth and satisfies ([?/)) elsewhere. In that case, S, as determined globally by 
((28)) . will still remain C^) smooth, by elliptic regularity. 

A natural weak formulation of l[27|) would be to look for p e satisfying 

L 

cl)^^K{p) + cl),x{p)pS^ dx = 0, Vc/. G Co°°(0, L), (34) 

where K[p) is a primitive for D{p). 

If we suppose that p has discontinuities at a; = s^, i = 1, ...,n, then the jump conditions for a 
solution of ([34]) satisfying the Neumann conditions at a; = 0, L are calculated to be 

{x{p)pS.,~D{p)p,){sf)^0 and K{p{sT)) = K(M)). (35) 

The first of (I35p implies, by the continuity of Sx, that the density gradients on either side of a 
discontinuity are coupled via 

D{p) _ D{p) + 

Candidate weak solutions can be constructed by solving, on each interval (s^, s^+i), the dynamical 
system 

G{pU^G{p)-p + Q, (37) 

for suitable Ci e M (see below), such that either p < p** or p > in each phase, and patching together 
the orbits to get a discontinuous density p on [0, L]. Once this is done, S is globally determined by 
solving the Neumann problem for (j28p . as already mentioned. 

For a pair (p, S) determined in this way, we can recover equation (j27p on each interval, as in 
Section 2.3.1. Thus, substituting ([28]) into (jST]) for each i, we arrive at 

n 

E 

1=1 



{||G(p) - 5 + C,||l,,((,^,,,^^)) - [(G - 5 + a)(G - 5 + G,).].} = 0, (38) 



where [-ji denotes the leap at a; = s^. 

It follows that G(p) = S* — Gi in each phase if, for example, the endpoints of adjacent orbits are 
chosen so that both the flux condition ([36]) and the jump condition [(G(p) -|-Gi)]i = hold. A global 
solution, p, constructed in this way is also a weak solution in the sense of (|34|) iff neighbouring Ci 
are chosen so that [K{p)]i — Vi. 

Remark. Steady states of the Stefan-problem formulation for ([3])- ([61), to be introduced in Section 5, 
are just special cases of such weak solutions for which the endpoints p{sf) take on a particular pair 
of a-dependent values; see Section 5.1 for details. 

3 Existence and uniqueness results 
3.1 Global existence for a < | 

In the case a = 0, a global-existence theorem for ([3]), ^ follows directly from the ideas of [B]. We 
now extend this result to cover all a G [Oi |)- 
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Theorem 3 Given smooth initial data (pq.Sq) satisfying < po < 1,0 < So < 1, the system (0], 
lO^; ^'^•s unique smooth global solution, (p, S), satisfying 0<p<l, 0<5<1, provided 
a < |. 

Proof. Since smooth data satisfying the Neumann boundary condition at a; = 0, L can be 
reflected about x = to give data on S^, and since ((21)-® are invariant under the transformation 
X I— > —X, it is enough to prove existence and uniqueness of smooth solutions given data on S^. 

For this, first note that (O can be used to write 

S={1- A)-V, (39) 

and that (1 — A)~^ is a bounded operator from L'^{S^) to H'^{S^). This can now be substituted into 
([3]), thus reducing our problem to the single nonlocal diffusion equation 

This equation can be solved using moUifiers, as in Section 15.7 of [H]- 

Specifically, for i?" data, we hit the right-hand side of (HDl) with a Friedrichs moUifier J^, replace 
p with Jep wherever it appears, and replace the initial data po with J^po- This results in an ODE in 
iJ", for Jep, such that the right-hand side is Lipschitz continuous in this quantity. Local existence 
and uniqueness of solutions follows by Picard's Theorem, as applied to Banach spaces [8] . 

As a consequence of the boundedness of (1 — A)~^, there is no obstacle to proceeding as in 
Chapter 15 of [9] to get the key estimate 

^^wpmIh^ < ciiiwmc^KWpemi. + 1). (41) 

With the aid of Sobolcv imbedding for s > |, this gives, as in Lemma 7.1, Theorems 7.2, 7.4, 
Chapter 15 of [5], a sequence p,^ ^ p e C{[0,T],C^{S^))nC°°{{0,T) x S^), such that p solves ^ 
on some time interval [0,T]. 

Next, Theorems 8.3, 9.6 and 9.10 of (Taylor III, Ch.l5), along with the Maximum Principle, 
imply that a continuation criterion for (|40p is that the chemotactic term be bounded in LP, some 
p > §, on finite time intervals, n being the space dimension. Thus, global existence of solutions will 
be established if we can show that the quantity 



^((l-.)(l-«p).f 



(42) 



is bounded on each [0,r]. For this, it is sufficient to show that ||ff is bounded on each [0,T], 
since Sx and Sxx are a priori uniformly bounded for all t, by the Maximum Principle applied to 

We proceed as in Lemma 2.2 of [3], and define an approximation of the sign function by 

as{z) = <j{z/6), for < 5 < 1, (43) 

with a a smooth and increasing function such that cr(0) = and 

cr(z) = sgn(z) for \z\ > zq, (44) 

some Zq > 0. 
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Then, setting abss{z) :— Jq crsiO we get abs5(z) — » \z\ as (5 — > 0, uniformly in z. Also note 
that 

, z > zqS , . 

0(i) , z<zo5. 



<^si^) = { 

Next, we differentiate ([3]) w.r.t. x, and mutiply by as{dxp)), which, upon integration, leads to 



It 



/ ahss{d^p) dx ^ Qs + 
Jo Jo 



as{d.,p)dUD{p)p.) dx, (46) 



where Qs comes from the chemotaxis term, and can be treated as in resulting in 

limlQal <Ci + C2 / \d^p\ dx. (47) 

For the remaining term on the rhs of (|46p . we get, writing v = p^, and integrating by parts, 

!^cT',iv)i-Dip)vl-D'ip)v^v')dx 

< a's{v)i-Dip)vl + \D'{p)\\v,\v^) dx 

(48) 

< C a'siv)v^ dx 

< CS\ 

where we have used the condition crsiO) = to kill the boundary terms, along with a'^ > 0, equations 
(|45)) . the elementary inequality AB < CqA^ + B^/ACq, and the fact that D{p) > 1 - |a > 0. 
Thus, taking the limit (5 — > in we see that 

^ \d,p\ dx<Ci + C2 \d^p\ dx. (49) 

An application of Gronwall's inequality now shows that the continuation criterion is satisfied. 
Finally, uniqueness can be proved with the aid of an analogue of ^1)) for the difference of two 
solutions, along with a Gronwall estimate for the L^-norm of p □ 

One might also expect global existence of smooth solutions for a > |, provided the initial den- 
sity profile were uniformly outside (either above or below) and xo were small enough. This can 
actually be proved for data close enough to a uniform steady state, as we now demonstrate. 



3.2 Global existence for a > |, xo small, and small data 

It turns out that we can use the comparison argument of Theorem 1 to obtain a global-existence 
theorem for ([3])- ([6]) when a > |, provided xo is small enough and the initial density profile, po{x), 
is sufficiently far from the unstable interval Jq. We can obtain results for both of the cases po{x) < 
p^ Mx G [0,£] and po{x) > p^ \/x E [0,i], but for clarity we will simply concentrate on the case 
Po{x) < p^ in what follows. 

First note that a local-in-time solution is guaranteed by the J^-method used previously, and that 
for global existence we merely need to prevent the solution from hitting p^ . 

Thus, let p be given, and pick pg-^ > p such that D[p) > Si for p < ps^. Next, pick an initial 
datum poix) such that avg(/Oo) = P, and such that maxpo(2:) < PSn and introduce a smooth. 
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modified diffusivity D*{p) which is equal to D{p) for p < ps^, and which is greater than for 

P> PSi- 

Equation ([3])* is then defined by replacing D{p) with D*{p) in the right-hand side of ([3]). Global 
existence of a solution {p, S) to ([3])*, (O, ([S]) follows by previous arguments, and it remains to show 
that p < ps^ yt, provided xo and ||po — p\\oc are chosen small enough. 

For this, let 

^2 = Xo max (1 - - ap)p, (50) 

0<p<p'> 

and define 

e = - 62. (51) 

If Xo is chosen so small that e > 0, then, in the same way as for a < |, we obtain the L'^-decay 
estimate 

\\p-pUt)<\\po-p\\2e-^'/\ (52) 
along with a sequence C" S [j^i + 1], such that 

\\p~p\UC)<C{Si,xo)\\po-ph (53) 

An i/^-estimate of the form ([T9l) again holds, whereby we emphasise that the right-hand side is 
Odlpo^plb), and it therefore follows from ([53|) and the comparison argument of Theorem 1 that 

\\p~p\Ut)<CiSi,xo)\\po-ph. (54) 

Thus, p < psi yt if \\po — p\\oo is chosen small enough, and consequently {p,S{p)) solves ©-([HI) for 
all time. 

Given a globally existing solution, convergence to the uniform steady state can be also be proved 
by the method of Theorem 1, provided xo is sufficiently small. 

The case p > p'* is handled analogously, and we therefore have, in summary 

Theorem 4 Given a smooth initial datum po{x) satisfying either po{x) < p'' Vx G [0, i] or po{x) > 
p^ Vx G [0, L], and letting p :— avg{po), equations (0)-® have a unique, local-in-time solution 
{p, S), which continues to exist globally if both ||po ~ p||oo cind Xo o'^e small. The smallness required 
of Xo depends on p and ||po ~ pIIoo- Furthermore, given any globally existing smooth solution which 
misses la, long-time exponential L°° -convergence to the uniform steady state holds, provided xo is 
sufficiently small. 

It should be noted that this result does not rule out the possibility of a globally existing solution 
with a > I approaching some non-uniform steady state outside as t ^ 00; in this regard, see the 
simulations in the next section. 

4 Numerics for the discrete model 

In this section we numerically solve the Neumann problem for ([3])- ([6]) in the high-adhesion regime 
by means of an explicit finite-difference scheme on a uniform spatial grid of n mesh points. 

The discretisation of the diffusion term in ([3]) is obtained by setting xo = in the right-hand side 
of ll]), while the chemotactic term is discretised by means of a simple, 0(/i)-accurate, upwinding 
scheme. Thus, in order to obtain numerical stability in well-posed regions, we have chosen a method 
of lines which is slightly different from the original discrete model ([ij. To complete the scheme, 
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we solve the elliptic equation ^ at each time step via the usual discrete Laplacian, together with 
Matlab matrix inversion. 

In what follows, we are particularly interested in observing how our numerical solutions change 
in the vicinity of as n increases, since this is essentially the same as asking in what sense ^ is 
the continmmr limit of ([1]) when a > |. 

4.1 Singular aggregation patterns 

We begin by choosing a — 0.95, xo — 16, and discretising a small initial density profile on a domain 
of length L = 8, using a grid of n = 400 spatial points. Evolving this data with our numerical 
scheme results in the sequence of snapshots displayed in Figure 2. 

Clearly, the effect of chemotaxis is to draw the solution towards , and once it has penetrated 
sufficiently far into the unstable region, a small number of fine oscillations quickly develop. Sub- 
sequently, mass is sucked into the central oscillatory region via a combination of chemotaxis and 
backward diffusion, while positive diffusion flattens out the density profile on either side. Eventu- 
ally, after a slow process of coarsening in which the fine oscillations disappear, we are left with a 
single, sharp-edged plateau, which presumably represents a steady-state weak solution of ([ll)-®, as 
constructed in Section 2.3.2. 

One further point to note here is that the values between which p jumps at the plateau edge are 
very close to those observed in [2] for the case xo = 0, where we also saw oscillations, as well as plateau 
formation through coarsening. In that paper the fact that the jump values appeared to depend only 
on a, and not on the initial data, was explained by the existence of a unique heteroclinic cycle for 
an 0{h^) modified equation derived from ((!]). We claim that the same argument goes through in the 
case Xo > 0, since the chemotactic terms produce a higher-order correction due to the smoothing 
properties of ([6]). The relevant saddle-point values for a = 0.95 are, in notation that will be used 
again in Section 5, (pi,p2) = (0.055,0.99) (see Figure 8 of [2]). 

Next, we repeat the simulation of Figure 2, but this time with n = 800 spatial points; the 
resulting snapshots are depicted in Figure 3. Notable differences with respect to Figure 2 are that 
the oscillations appear a little earlier, and at a slightly lower density level, and that the jumps levels 
at large times appear even closer to pi and p2- 

Finally, we repeat the simulation once more using n = 1200 spatial points; the solution is plotted 
in Figure 4. Note that the level at which oscillations arise is now even earlier, and closer to the 
lower end of /„, but that the rate of convergence to (presumably) p'^ is exceedingly slow w.r.t n. 
The jump levels have closed in even further on pi and p2- 

4.2 Macroscopic coarsening 

Putting aside microscopic oscillations (which are related to ill-posedness) , for the moment, a phe- 
nomenon exhibited by other, well-posed chemotaxis models is that of macroscopic coarsening, 
whereby a large aggregation region attracts a smaller one (see, e.g., jZl |4j). Such behaviour is 
in fact also exhibited by our model ©-([I]), as evidenced by the simulation of Figure 5, in which 
a wide plateau region absorbs a much narrower neighbour, resulting in a (quasi?) steady state at 
large times. 

This brings us to another phenomenon associated with chemotaxis equations subject to the 
Neumann condition, which is that a single, asymmetrical plateau will tend to move (perhaps very 
slowly) towards the boundary as < ^ oo. Unfortunately, our numerical code is not accurate enough 
to say definitively which way (if any) the plateau in Figure 5d is moving; one might hazard a guess 
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Figure 2: Evolution of a small initial density profile, using n = 400 spatial points. Data shown at 
(a) t = 0, (b) t = 1.3, (c) t = 1.38, (d) t = 1.8, (e) t = 6, (f) t = 7. The parameter values are 
a — 0.95, xo = 16, i = 8, and the boundaries of the unstable region are marked with dotted 
lines. 
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Figure 3: Evolution of the same data as in Figure 2 using n = 800 spatial points. Data shown at 
(a) t = 0, (b) i = 1.1675, (c) t = 1.255, (d) t = 1.805, (e) t = 7.2, (f) t = 7.6. The parameter values 
are a = 0.95, xo = 16, i = 8. 
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Figure 4: Evolution of the same data as in Figure 2, using n = 1200 spatial points. Data shown at 
(a) t = 0, (b) i = 1.0911, (c) t = 1.1356, (d) < = 2, (e) < = 8.6778, (f) t = 9.344. The parameter 
values are a ~ 0.95, xa = 16, L = 8. 
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Figure 5: Evolution of an initial density profile having one large and one small plateau. Data shown 
at (a) t = 0, (h) t = 8, (c) t = 13, (d) t = 18. The parameter values are a = 0.95, xo = 16, L = 8. 



that the plateau remains where it is, due to the fact that there is essentially a Dirichlet condition 
on either side of the jump locations. 

4.3 Smooth, non-uniform steady states 

As mentioned below the statement of Theorem 4, our analytical results allow for the possibility that 
a solution of ([3|)-(l6|) with a > | might approach a smooth non- uniform steady state, avoiding /„, as 
t — > (X). Some numerical evidence for this is presented in Figure 6, which depicts overlayed snapshots 
of a low-mass, high-adhesion solution converging to a bell-shaped steady state. 
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Figure 6: Evolution of a small initial density profile towards a non- uniform steady state below 
la = [0.361,0.973]. Data shown at t = 0,2,4,6,8,10 and 12, such that the central maximum 
increases with time. The parameter values are a = 0.95, xo = 16, L = 8. 
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4.4 The question of critical mass 

For some chemotaxis models, such as the Keller-Segel equations in R^, to take a weh-known example, 
there is a bifurcation phenomenon, such that solutions exist for all time (and also disperse) if the 
mass is below some critical value, while blow up (formation of Dirac deltas) occurs in finite time 
otherwise [5]. We do not, however, expect our adhesion/chemotaxis model to exhibit quite this kind 
of bifurcation, since backward diffusion and volume filling have the effect of stabilising even very 
slender aggregations. 

A numerical example of this is shown in Figure 7. Here, we took an initial datum with a thin 
high-density region, and a mass so small that the stability condition ((26)l is satisfied. Despite the 
uniform steady state being locally stable, the aggregation region appears to persist, such that the 
overlayed snapshots are visually indistinguishable. We would also expect to see analogous behaviour 
in an appropriate 2-d version of ([SI)- ([6]) in the high-adhesion regime. 

A slightly different, but related, question one can ask is whether a low-mass aggregation satisfying 
(|26p can be obtained by evolving an initial datum which lies below /„, but has, say, a narrow spike 
almost touching . Despite numerous attempts, we have been unable to achieve this numerically; 
the spike always collapses almost immediately. Thus, there is evidence that (|26|) implies stability 
with respect to perturbations which remain below . 

5 Stefan problems 
5.1 Formulation 

In Section 4.1 we touched on an observation made in [2 for the special case xo = 0; namely that 
large-time plateau values in solutions of ([T]) seem to be essentially unique, for a given value of a, 
and we noted that such uniqueness is also expected to hold for xo > 0, as a consequence of elliptic 
regularity. The observation of [2] subsequently led to the idea that a Stefan-problem framework, in 
which solutions are allowed to jump between unique plateau values pi{a) and p2{oi) (as calculated 
in [2]), but are elsewhere smooth, might be an appropriate way of treating fS]) as the limit of ^ in 
the high-adhesion regime, and it did indeed prove possible to develop an (at least partial) existence- 
and-uniqueness theory for such problems [T]. Continuing in this vein, we will now proceed to write 
down a Stefan-problem formulation for (O-® in the simplest possible case. 

Imagine, then, that we are given a small initial density profile po below the unstable region (i.e., 
such that Pq{x) < p^ Vx), and imagine that we evolve this data via ([3])-(l6|) until p hits p'' at some 
point Xc and time tc- The idea now is to continue the solution past tc by means of a singular, three- 
phase Stefan problem, whereby we introduce a high-density middle phase at Xc, which is initially 
and instantaneously a zero- width spike jumping up from pi to p2 (and back down again), and which 
will subsequently fatten up as mass is drawn in from the low-density left- and right-hand phases. 
In each of the left and right phases we impose the Dirichlet condition p = pi at the boundary with 
the middle phase, and in the middle phase we demand that p — p2 ai the left and right boundaries, 
which will be denoted by si{t) and Sr{t). The density in each phase evolves according to the 
moving boundaries si{t) and Sr{t) evolve according to appropriate Rankine-Hugoniot conditions, and 
finally, since all of this is rather difficult to explain in words, we refer the reader to the simulation 
of Figure 8 (discussed below) for further clarification. 

To be mathematically explicit, for t > tc we wish to solve in each phase the adhesion/chemotaxis 
equation 




(55) 
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Figure 7: Overlayed snapshots of a low-mass aggregation at essentially steady state. Data shown at 
t = 0, 0.2222, 0.6667, 1.1111. The parameter values are a = 0.95, xo = 16, L = 8. 
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subject to the boundary conditions 



dp 
dx 



— = at X — Q and p = pi at x = (t) (56) 



for the left-hand phase, 

p = P2 at x — sl{t) and x = {t) (57) 
for the middle (high-density) phase, and 



p — Pi at X = s^{t) and ^ = at x = L (58) 

ox 

for the right-hand phase. 

The chemoattractant concentration S is obtained by globally solving the Neumann problem for 

AS^S-p., (59) 

which is well-posed despite p having two jump discontinuities, by elliptic regularity, as noted earlier. 
The moving boundaries si{t) and Sr{t) are governed by the pair of Rankine-Hugoniot conditions 

^ = ~i{Dipi)px - x{pi)PiSx){sY) - {D{p2)Px - x{p2)p2Sx){st))/{pi - P2), (60) 

dSr 

-{{D{p2)px - x{P2)p2Sx){s';) - [D{pi)px - x{pi)PiSx){st))l [P2 ^ Pi), (61) 

which guarantee local conservation of mass. 

Unfortunately, equations ([55|) - (pT|) have proved to be analytically intractable when subject to the 
singular initial condition si{tc) = Sr{tc)- In particular, we have been unable to prove the (plausible) 
conjecture that the Dirichlet condition p = pi at Sj~ and a priori holds the density below p^ in 
each of the low-density phases for some short time. 

Nevertheless, we can at least attempt to solve these equations numerically, as we do below, 
whereby our attention will be focused on three questions: 

1. Do sensible- looking solutions exist? 

2. As a matter of principle, are solutions close (in some weak sense) to the oscillatory solutions 
obtained by discretising ©-(El) directly? 

3. Might it be more computationally efficient to solve the Stefan problem than to discretise ([3])- ([51) 
directly? 

5.2 The rescaled model; numerical solutions 

Following the approach of [1] , we solve the chemotaxis equation in each given phase by rescaling the 
spatial coordinate so as to fix the relevant moving boundary (or boundaries). 
In the left-hand phase this results in 

dp 1 d ( dp\ sidp 1 d ( dS\ 
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for X e [0, 1], in the middle phase we get 

dp 1 ^ ( r\i \'^P\ [srX + [1 — x)si) dp 1 d ( dS 



dt {sr - sif dx y^'^P^Tx ) + (..-.0 Tx " T^—WTx y'^p^p^ ) ' (^^^ 

for X e [0, 1], and in the right-hand phase 

dp 1 d ( ^dp\ , {{i-x)sr) dp 1 d ( ds\ 



dt {L — s,.)2 9a:: \ dx J {L — s,.) dx (L — s,.)^ 9x \ 9a; 
for X e [0, 1]. 

The rescaled Rankine-Hugoniot conditions take the form 

^ = - f- iDipi)p, ~ x{pi)piS.) {!-) - , ^ . {D{p2)p. - X{P2)P2S,) (0+)) /(pi - p2), 

(65) 

^ - - i i ^ ^ (^(P2)P. - X(P2)P25,)(1-) - ^ (i^(pi)p. - X{Pl)PlS.) (0+)) /(P2-Pl). 

(66) 

Each of equations is solved on a uniform grid in essentially the same way as in Section 

4, while (|65 p - (|66|) are solved (explicitly) using one-sided, second-order-accurate finite differences for 
the gradients. Since the three phases are generally of different lengths, this entails that the numerical 
approximation of p lives on a globally non-uniform grid. In order that we can nevertheless solve (|59p 
using the discrete Laplacian, p is interpolated onto a globally uniform grid at each time step. 

To obtain the simulation shown in Figure 8, we used the numerical method of Section 4 to evolve 
the initial data of Figure 4 (with n = 1200 spatial points) until the solution hit p = p^ at tc — 0.8325, 
and then continued the solution via the Stefan-problem algorithm just described, such that there 
arc 100 spatial points in each of the three phases. 

Note that the solution is nice and smooth away from the moving boundaries, and that that 
Dirichlet conditions at si and Sr hold the solution below p^ for all time in the low-density phases. 
Also, since the gradient at sj' and is large just after ic, the middle phase gains mass very quickly 
for a short time; subsequent to this, there is a slow approach to the kind of weak steady-state solution 
of ©-dSl) seen in Figures 2-4. 

In Figure 9 we overlay the simulations of Figures 4 and 8 in order to compare the Stefan-problem 
approach with that of direct discretisation. We see that, away from the central oscillatory region, 
there is always good agreement between the solutions, but that towards the middle of the domain 
there is a significant discrepancy shortly after tc, due to the fact that, with the direct approach, the 
density has to push a considerable distance into la before oscillations set in (even when n — 1200), 
thus creating a short delay. Also, the solution of Figure 4 unfortunately gains a little mass during 
the course of the simulation. Despite this. Figures 2-4 and 8, taken together, do seem to indicate 
that the Stefan problem is the correct (weak) limit of ([T]). 

One advantage of the Stefan-problem framework is that the microscopic oscillations seen in the 
direct approach are avoided, and one can get convincing solution behaviour using relatively few 
spatial points. However, one disadvantage of the rather obvious numerical method we used (and 
which we certainly don't claim to be the best) is that, due to the spatial rescalings and the large 
initial gradients at sJ' and , the parabolic and hyperbolic CFL conditions demand a very short 
time step until the high-density phase has attained a considerable thickness. 
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Finally, one can of course imagine solutions of the three-phase Stefan problem in which chemo- 
taxis is so strong that the solution in one of the low-density phases rises up to hit once again. In 
that case, another spike can be inserted at the point of contact, and the solution continued via the 
appropriate five-phase Stefan problem, and so on, ad infinitum. 
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(a) 



(b) 




Figure 8: Here we show what happens when the n — 1200 solution of Figure 4 is continued via the 
Stefan problem after hitting at t = 0.8325 (see (a)). In each phase there is a uniform mesh of 
100 points. Snapshots (b)-(f) are taken at the same times as in Figure 4, namely (b) t = 1.0911, (c) 
t ^ 1.1356, (d) t = 2,{e)t = 8.6778, (f) t = 9.344. 
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Figure 9: Initial data for the Stefan problem of Figure 7 (a), together with synchronised, overlayed 
snapshots from Figures 4 and 7 (b)-(f). 
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